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RESEARCH  OVERVIEW 


The  work  under  this  contract  is  carried  out  in  four  areas.  These  are  con¬ 
cerned  with  CAF  system  structure,  process  simulation  and  modeling,  equipment 
modeling  and  simulation,  and  scheduling. 

Our  CAF  system,  named  CAFE  (Computer-Aided  Fabrication  Environment)  is  in  con¬ 
stant  use  in  the  MIT  IC  Laboratory.  It  currently  runs  on  a  VAX  785,  but  will 
shortly  be  ported  to  a  Sun  Microsystems  computer.  At  this  time  the  data  base 
will  be  converted  so  it  uses  INGRES  running  on  the  Sun. 

A  data  model  for  VLSI  fabrication  facilities  is  currently  under  development. 
There  are  various  type  extensions,  and  the  schema  captures  several  important 
aspects  of  pJLant  and  process  management. 


Progress  has  been  made  on  the  definition  of  a  language  for  specifying  process 
flows.  This  work,  done  in  loose  cooperation  with  the  Berkeley  CIM  effort,  is 
no.  yet  at  a  stage  where  it  can  be  used.  However,  to  facilitate  operation  in 
the  IC  Laboratory  a  very  abbreviated  and  preliminary  version  of  a  flow  lan¬ 
guage  h^s  been  implemented,  in  the  sense  that  several  interpreters  have  been 
written  that  act  on  the  CMOS  baseline  process  under  development.  The  language 
effort  is  based  on  the  two-stage  generic  process-step  model  which  has  been 
disptissed  in  previous  reports.  This  model  also  forms  the  foundation  for  the 
flow  language  effort  at  Berkeley  and  Stanford,  as  reported  at  the  informal 
Workshop  on  Process  Specification  held  at  the  conference  on  Advanced  Research 
in  VLSI,  March  24,  1987,  at  Stanford  University. 

Two-dimensional  localized  thermal  oxidation  has  been  modeled  using  a  visco¬ 
elastic  boundary  element  method.  This  model  predicts  higher  than  expected 
stress  during  oxidation,  and  hence  possibly  plastic  deformation. 

There  is  a  need  for  an  interchange  format  for  wafer  state..  We  have  developed 
a  suite  of  routines  to  work  with  the  propose^,  s-tamterff  profile  interchange 
format. 

-‘The  equipment  modeling  task  is  a  new  one.  The  objective  is  to  develop 
executable  machine  models  to  be  used  for  machine  simulation.  The  models  will 
deal  with  both  nominal  values  and  variances.  A  matrix  experiment  will  be 
devised  to  test  such  models  for  LPCVD. 


A  hierarchical  framework  for  scheduling  has  been  devised,  in  which  different 
levels  deal  with  different  time  horizons.  .  The  lowest  level  is  the  actual, 
detailed  set  of  schedules.  An  equipment  fih&gervation  system  for  this  lowest 
level  has  been  written  and  incorporated  into  CAFE,  and  is  now  in  daily  use. 
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CAF  SYSTEM  STRUCTURE 


Last  year  an  initial,  rudimentary  CAF  system  was  designed  and  implemented. 

This  included  a  personal  notebook  and  data  structures  for  the  upcoming  fabri¬ 
cation  equipment  installation. 

Our  present  computer  hardware  consists  of  a  VAX  785  with  16  Mbytes  of  main 
memory,  five  450-Mbyte  disks,  GCR  tape,  1600  cpi  tape,  two  laser  printers, 
four  phone  lines,  and  an  Ethernet  port  to  the  MIT  network.  Additionally  there 
are  eight  11/73-based  terminal  concentrators,  six  with  48  RS232  ports,  one 
with  16  ports  and  one  spare  for  maintenance.  RS232  cables  have  been  installed 
throughout  our  building,  and  numerous  user  terminals  and  PCs  have  been  con¬ 
nected. 

Because  of  our  success  in  getting  faculty,  staff,  and  students  to  use  our 
computing  facilities,  and  because  of  the  increased  use  of  our  CAF  system  in 
actually  operating  the  IC  Laboratory,  our  VAX  785  (CAF.MIT.EDU)  is  heavily 
overloaded.  Ve  must  augment  our  facilities.  We  evaluated  several  possibili¬ 
ties  and  have  initiated  the  purchase  of  two  Sun  Microsystems  3/280  computers. 
One  of  these  will  be  used  for  the  actual  running  of  the  CAF  system  CAFE  for 
the  IC  Laboratory  and  the  other  will  be  reserved  for  development  of  CAF  system 
software.  We  also  have  ordered  RTI  INGRES  for  use  on  the  Sun  computers. 
Delivery  of  this  equipment  is  expected  during  the  first  half  of  1987. 

We  have  made  substantial  progress  in  the  development  of  a  data  model  and 
schema.  Our  data  architecture  being  developed  is  similar  to  the  Multibase 
system  developed  at  CCA.  This  provides  a  uniform  query  interface  to  data 
residing  in  multiple  autonomous,  heterogeneous  data  bases.  Our  current  data 
base  is  distributed  across  two  relational  systems,  university  INGRES  and 
PRELUDE  and  a  hierarchical  file  system.  PRELUDE  is  a  fast,  lightweight  UNIX- 
style  data  manager.  Our  system  is  very  modular  and,  thus  far,  it  has  been 
easy  to  incorporate  new  DBMSs  into  the  system  as  well  as  move  data  from  one 
data  manager  to  another. 

Our  data  model  is  the  functional  model  with  support  of  extended  data  types 
including  various  temporal  types  as  well  as  inexact,  interval,  and  null 
values.  The  schema  captures  several  important  aspects  of  plant  and  process 
management;  fabrication  facilities  and  equipment,  users,  equipment  reserva¬ 
tions,  iots,  lot  tracking,  wafers,  process  flow  descriptions,  wip  tracking, 
and  lab  activity  information. 

We  have  developed  a  generalized  forms  based  user  interface  program,  called 
fabform.  This  single  program,  when  called  with  a  parameter  file,  produces  a 
terminal  display  and  allows  a  user  to  move  from  field  to  field  and  enter  data. 
The  type  and  content  of  the  screen  display  is  specified  by  an  ASCII  file  which 
is  referenced  by  data  included  in  the  parameter  file.  The  form  may  have  arbi¬ 
trary  length  and  the  user  can  scroll  up  or  down.  At  present,  user  interface 
commands  are  much  like  EMACS  commands.  When  the  user  exits,  or  ”Ves  the 
data,  an  updated  parameter  file  is  written. 
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MODULAR  PROCESS 


The  implementation  of  the  first  version  of  MASTIF  (MIT  Analysis  and  Synthesis 
Tool  for  IC  Fabrication)  was  completed  in  May  1986.  A  menu-  and  window- 
oriented  program  has  been  developed  as  the  first  step  in  meeting  the  need  for 
an  integrated  process  design  system.  MASTIF  includes  aids  for  process  speci¬ 
fication  and  simple  process  verification,  and  provides  interfaces  (both  front 
and  back  ends)  to  the  SUPREM-III  process  simulatoL  and  the  MINIMOS  device  sim¬ 
ulator.  Drawing  on  analogies  from  other  areas  of  VLSI  design,  a  blueprint  for 
future  development  of  tools  beyond  simple  simulation  has  been  completed  [1], 
MASTIF  was  successfully  ported  to  run  on  a  true  workstation,  the  Vaxstation-II 
under  VMS.  The  program  is  now  in  its  second  phase  of  development  and  in  the 
process  of  being  ported  to  a  UNIX  environment. 

We  have  continued  our  effort  on  modeling  two-dimensional  (localized)  thermal 
oxidation.  Thermal  oxidation  of  silicon  involves  the  diffusion  of  oxidant 
species  from  the  gas-oxide  interface  to  the  oxide-silicon  interface,  and  the 
transport  of  newly  formed  oxide  away  from  the  latter.  Under  suitable  formula¬ 
tions,  it  can  be  shown  that  the  diffusion  process  is  a  Laplace  problem  and  the 
viscoelastic  flow  of  oxide  is  a  biharmonic  problem.  For  these  boundary  value 
problems,  the  unknown  boundary  parameters  can  be  obtained  from  the  known 
boundary  conditions  without  calculating  the  interior  solutions.  The  diffusion 
problem  is  solved  with  a  standard  boundary  element  method  (BEM)  for  potential 
problems.  A  generalized  viscoelastic  BEM  has  been  developed  to  model  the 
oxide  flow.  Utilizing  constant-velocity  kernel  functions,  this  viscoelastic 
BEM  can  deal  with  a  wide  range  of  stress  relaxation  times,  covering  elasto- 
static  deformation  and  incompressible  creeping  flow.  Our  approach  achieves 
simplicity  and  efficiency  by  solving  a  two-dimensional  problem  as  line  inte¬ 
grals  on  the  boundaries.  Simulations  of  Local  Oxidation  of  Silicon  (LOCOS) 
structures  indicated  that  stress  created  during  oxidation  could  be  extremely 
high,  particularly  when  the  mechanical  barrier  effect  of  silicon  nitride  mask 
is  included.  This  suggests  that  both  silicon  nitride  and  oxide  flow  or  plas- 
ticly  deform  more  readily  than  assumed.  Stress-induced  retardation  of  reac¬ 
tion  rate  and  diffusivity  of  oxidants  are  also  studied.  In  addition  to  an 
overall  lowering  of  stress,  the  shape  of  the  oxide  changes  significantly  when 
such  nonlinear  effects  were  included. 

One  important  finding  of  the  MASTIF  project  has  been  to  note  the  need  for  a 
uniform  representation  of  wafer  and  device  structure  information,  both  for  use 
by  individual  tools  in  a  complete  design  system,  and  for  interchange  between 
different  simulation  sites.  Development  of  general  analysis  tools  and  inter¬ 
faces  requires  a  central,  agreed-upon  representation  of  the  structures  these 
simulators  and  associated  tools  manipulate.  We  have  been  heavily  involved 
in  standards  work  related  to  a  Profile  Interchange  Format  (PIF) ,  and  are 
implementing  a  library  of  routines  for  accessing  these  structures  based  on  the 
PIF  (2).  Specific  tools  and  components  urder  development  include  SNC,  a  local 
"database”  form  for  PIF  storage,  PIFLIB,  a  library  of  routines  for  tool  inter- 
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EQUIPMENT  MODELING 


This  project  is  a  new  one,  not  included  in  previous  reports. 

VLSI  machine  modeling  has  taken  a  concrete  form  during  the  past  several 
months.  The  thrust  of  the  work  is  to  combine  analytical  modeling  with  matrix 
experimental  approaches  and  to  provide  an  executable  program  which  will 
facilitate  the  following: 

1.  Off  line  quality  control  to  determine  the  point  in  operating  space  which 
provides  the  greatest  robustness  against  variations  in  process  parameters  and 
therefore  yields  the  most  consistent  results. 

2.  On  line  quality  control  used  to  tune  a  process  based  either  on  measure¬ 
ments  made  after  a  previous  run  or  on  in-situ  measurements. 

Matrix  experimentation  is  a  collection  of  methods  wherein  some  or  all  of  the 
relevant  process  variables  are  varied  simultaneously  in  an  experiment  and 
information  is  extracted  from  the  results  statistically  [1,2].  In  contrast  to 
conventional  single-variable  experimentation,  the  matrix  approach  offers  a 
tremendous  economy  of  experimental  effort.  This  is  crucial  in  a  production 
environment  as  the  interruption  to  work  flow  must  be  kept  to  a  minimum,  while 
in  a  research  environment  it  is  useful  in  order  to  optimize  a  new  process  as 
quickly  as  possible.  Matrix  experimental  techniques  have  been  employed  with 
great  success  on  VLSI  processes  at  AT&T  Bell  Laboratories  for  approximately 
the  last  six  years  [3,4], 

Our  process  modeling  effort  will  begin  by  utilizing  analytical  and  experience 
based  background  to  define  the  process  variables  and  their  values  for  a  matrix 
experimental  approach.  The  analytical  work  will  be  based  both  on  process 
models  available  in  the  literature  and  on  simple  physical  models.  Future 
extensions  of  the  work  will  attempt  to  more  closely  couple  the  analytical  and 
experimental  approaches  by  using  matrix  experimentation  to  verify  analytical 
models,  specify  numerical  values  for  the  analytical  models,  and  also  to 
improve  analytical  models. 

The  first  process  selected  for  the  machine  modeling  effort  is  LPCVD  of  poly 
silicon.  Optimization  of  thickness  uniformity  across  a  wafer,  between  wafers 
in  a  single  lot,  and  between  runs  will  be  the  goal. 


References: 

[1]  G.  Box,  W.  Hunter,  and  J.  Hunter,  "Statistics  for  Experimenters  and 
Introduction  to  Design  Data  Analysis  and  Model  Building,"  Viley  1978. 

[2]  G.  Taguchi,  "Introduction  to  Quality  Engineering,"  Krauss  International 
Publications,  White  Plains,  NY,  1986. 
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SCHEDULING 


Research  during  this  period  focused  on  three  activities:  (1)  studying  the 
integrated  circuit  fabrication  process  at  a  systems  level,  (2)  formulating  a 
mathematical  model  of  an  integrated  circuit  fabrication  facility,  and  (3) 
developing  simulation  and  scheduling  software. 

The  effort  to  define  the  scheduling  problem  continues.  We  are  concentrating 
on  using  the  MIT  laboratories  as  case  studies.  Mathematical  and  simulation 
nodels,  described  below,  are  being  based  on  what  we  learn  here. 

As  a  mathematical  model  of  an  IC  fabrication  facility,  the  multiple  time  scale 
decomposition  under  development  shows  great  promise.  A  new  basic  model  is 
being  investigated  which  will  help  us  refine  and  better  justify  the  tentative 
mathematical  results  we  have  developed  thus  far  on  hier archied' scheduling. 

In  this  approach,  the  scheduling  algorithm  is  divided  into  a  set  of  levels 
which  correspond  to  classes  of  events  that  are  distinguished  by  their 
frequencies.  At  each  level,  two  kinds  of  calculations  are  performed:  small 
linear  programs,  to  determine  frequencies  of  higher  frequency  (lower  level) 
events;  and  simple  combinatorial  optimizations,  to  determine  exact  times  for 
the  events  of  that  level,  whose  frequencies  have  been  calculated  at  higher 
levels. 

Software  which  will  implement  a  multiple  time  scale  decomposition  approach  to 
hierarchical  scheduling  is  under  development.  A  simulation  is  also  under 
development.  The  scheduling  software  will  first  be  tested  with  the  simula¬ 
tion,  and  then  used  to  run  the  laboratory. 

An  electronic  machine  reservation  system  is  also  under  development.  In  its 
initial  version,  it  is  essentially  an  electronic  sign-up  sheet  for  equipment. 
It  will  later  be  used,  after  modifications,  as  the  lowest  level  of  the 
hierarchical  scheduler. 
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A  VISCOELA  jTIC  BEM  FOR  MODELING  OXIDATION 


Thye-Lai  Tung,  Jerome  Connor,  and  Dimitri  A.  Antoniadis 
Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts  02139 

Abstract 

A  viscoelastic  boundary  element  method  has  been  developed  to  model  the  motion 
of  silicon  dioxide  and  silicon  nitride  during  thermal  oxidation  of  silicon.  This  tech¬ 
nique  uses  Kelvin’s  solution  reformulated  according  to  the  correspondence  principle 
on  viscoelasticity.  Constant-velocity  loading  is  chosen  to  ensure  smooth  variations  in 
displacement  and  stress  behavior  for  a  wide  range  of  relaxation  times. 


1  Introduction 

A  major  problem  in  modeling  thermal  oxidation  of  silicon  is  the  moving  boundaries.  Driven 
by  the  conversion  of  silicon  to  silicon  dioxide  of  increased  volume,  boundaries  change  shape 
drastically  during  a  local  oxidation  of  silicon  (LOCOS)  process  step.  Conventional  methods 
such  as  the  finite  element  method  (FEM)  require  a  mesh  to  subdivide  the  simulation 
domain.  To  apply  such  methods  to  thermal  oxidation,  one  must  first  develop  computer 
codes  to  regenerate  the  mesh  automatically  and  optimally  as  the  oxide  changes  shape  at 
every  time  step  [1,2,3).  After  that,  one  may  have  to  deal  with  the  transfer  of  stress  history 
from  the  old  mesh  to  the  new  one,  depending  on  what  oxide  flow  model  is  used.  This  is 
another  difficult  problem  if  one  wishes  to  minimize  numerical  diffusion  of  stress  distribution 
due  to  regriding. 

’  In  contrast  to  the  FEM,  the  boundary  element  method  (BEM)  does  not  require  a 
mesh  in  general  since  it  does  all  calculations  on  the  boundary  but  not  in  the  interior.  Be¬ 
cause  segments  need  not  be  regenerated,  keeping  track  of  stress  history  on  the  boundary 
is  straight  forward.  However,  previous  efforts  on  applying  the  BEM  to  model  oxide  mo¬ 
tion  as  viscoelastic  flow  suffered  from  some  drawbacks.  In  Matsumoto’s  foi  mutation  [4], 
domain  calculations  were  needed.  Although  restrictions  on  the  grid  were  not  as  stringent 
as  those  in  the  FEM,  excessive  computation  time  was  evident.  Simulation  setup  might 
still  be  difficult  for  complex  structures.  Isomae  attacked  the  problem  differently  by  using 
a  Laplace  transform  technique  [5].  According  to  the  correspondence  principle  of  linear 
viscoelasticity,  we  may  solve  a  viscoelastic  problem  through  an  equivalent  elastostatic  sys¬ 
tem  in  the  Laplace  transform  space.  Boundary  element  techniques  for  elastostatics  do  not 
require  domain  calculations,  but  unfortunately  it  is  impossible  to  solve  Laplace  transforms 
analytically  in  any  practical  simulation  situations.  Numerical  Laplace  transform,  as  used 
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by  Isomae,  is  deemed  wasteful  because  a  new  solution  must  be  generated  at  every  time 
step. 

In  this  naper  we  describe  a  generalized  BEM  for  two-dimensional  linear  viscoelastic 
flow,  with  application  to  thermal  oxidation.  Instead  of  performing  the  Laplace  transform 
on  the  global  solution,  we  operate  it  on  the  fundamental  equations.  Our  kernel  functions 
are  derived  from  Kelvin’s  solution  used  in  elastostatic  BEM  applications  [6].  A  suitable 
excitation  function  is  added  to  produce  a  time-varying  body  force  that  causes  the  load 
point  to  move  with  a  constant  velocity  in  a  viscoelastic  medium.  This  approach  eliminates 
domain  calculations  and  numerical  Laplace  transforms.  Yet,  it  is  as  easy  to  use  as  those 
for  elastostatics.  The  main  disadvantage  is  that  boundary  conditions  cannot  be  satisfied 
for  all  times,  so  the  system  must  be  solved  periodically.  But  that  is  not  a  problem  for 
thermal  oxidation  because  a  new  solution  must  be  obtained  anyway  for  every  time  step. 
This  formulation  can  handle  all  possible  values  of  Poisson’s  ratio  and  a  wide  range  of 
stress  relaxation  times,  essentially  encompassing  viscous  incompressible  flow  and  elastic 
deformation. 

2  Numerical  Formulation 

We  implement  the  viscoelastic  BEM  using  the  so-called  indirect  formulation,  which  is 
also  known  as  the  classical  or  source  method.  Here  we  attempt  to  model  a  problem 
by  putting  “sources”  on  the  boundary  and  adjust  their  strength  so  that  the  fields  they 
generate  match  the  prescribed  boundary  conditions.  For  a  potential  problem,  the  sources 
are  simply  electrical  charges.  The  indirect  formulation  has  an  advantage  over  the  direct 
formulation  in  that  different  components  within  the  stress  tensor  or  displacement  vector 
can  be  obtained  more  readily. 

The  thermal  oxidation  process  actually  consists  of  two  tightly  coupled  processes,  namely 
oxidant  diffusion,  and  oxide  motion.  For  details  on  the  oxidant  diffusion  process,  boundary 
conditions,  and  numerical  implementation  of  the  integral  equations,  readers  are  referred 
to  our  previous  paper  on  the  same  subject  [7],  In  that  paper  we  use  an  incompressible 
viscous  flow  model  for  oxide  motion;  here  we  have  a  generalized  model  applicable  to  the 
oxide,  silicon  nitride  and  silicon  substrate. 

The  two  sets  of  integral  equations  for  viscoelasticity  are  shown  below: 


where  <r, 7- ,  u;,  and  p  are  the  stress  tensor,  displacement  vector,  and  soiure  density  respec¬ 
tively.  T  denotes  the  boundary.  o*-k  and  u*;  are  the  kernels;  their  actual  forms  are  given 
in  the  appendix.  For  compactness,  Einstein  notation  has  been  used.  To  get  the  surface 
traction,  we  apply  the  following  formula: 

Pi{x)  =  Oji(x)nj(x) 

where  rij  is  the  direction  cosine  of  the  surface. 

If  we  want  to  treat  silicon  nitride  as  an  elastic  material  with  stiffness,  we  have  to  modify 
the  boundary  condition  of  the  oxide-nitride  interface  from  a  free  surface  condition  p  =  0 
to  the  following: 

1.  displacement  vector  is  continuous  across  the  interface, 

2.  surface  tractions  of  the  two  materials  are  equal  and  opposite. 

To  ensure  stability,  the  equations  for  silicon  nitride  layer  and  oxide  region  must  be  solved 
simultaneously.  The  resulting  system  matrix  is  larger  but  banded. 

Likewise,  we  use  the  same  approach  if  we  desire  to  model  the  silicon  substrate  as  an 
elastic  foundation,  instead  of  treating  it  as  a  rigid  body.  Note  that  due  to  the  unique 
formulation  of  BEM,  boundary  conditions  at  y  =  — oo  need  not  be  specified  for  the  silicon 
substrate. 

In  viscoelasticity,  we  must  include  past  stress  history  in  the  present  time  step.  Consider 
p™,  the  normal  component  of  the  surface  traction  at  time  step  m.  The  two  different  relax¬ 
ation  terms,  as  defined  in  the  appendix,  are  kept  seperated  and  updated  in  the  following 
way: 

Pna  =  P„a  +p”_1eXp(-Ar/ra) 

Pn0  =  Pn0  +P™g'  exp{~Atm/T0) 

where  A tm  is  the  time  step  size,  ra  and  t0  the  characteristic  relaxation  time  constants. 

3  Simulation  Results 

We  will  demonstrate  how  stress  is  relieved  via  viscoelastic  flow,  *he  effect  of  the  silicon 
nitride  layer  on  the  oxide  shape  and  the  stress  distribution.  The  relaxation  time,  or 
equivalently  the  viscosity,  of  oxide  has  not  been  determired  accurately.  It  is  difficult  to 
do  so  because,  for  a  given  temperature,  the  viscosity  depends  on  the  quality  of  the  oxide, 
water  contents,  and  the  presence  of  other  impurities.  Most  of  the  recent  data  on  oxide 


viscosity  are  inferred  rather  than  measured  directly  in  experiments,  but  they  are  useful  as 
ballpark  figures.  In  any  case,  it  is  a  good  exercise  to  vary  the  relaxation  time  to  see  how 
it  affects  the  stress  distribution. 

Shown  in  Fig.  1  is  the  outline  of  a  semi-recessed  LOCOS  structure  plotted  for  every 
time  step.  In  this  simulation  window  of  I.6//m  wide,  the  silicon  nitride  mask  extends  from 
x  =  0  to  x  =  0.96/im  and  is  assumed  to  be  totally  flexible.  Initially  the  structure  has  a  pad 
oxide  thickness  of  200A.  Oxidation  is  carried  out  at  925°C  in  a  wet  ambient  for  3.4  hours 
to  get  a  final  field  oxide  thickness  of  5000A.  The  Young’s  modulus  and  Poisson’s  ratio  are 
taken  to  be  8  x  1011  dynes-cm-2  and  0.194.  Shown  in  Fig.  2a,  2b  and  2c  is  the  normal 
surface  traction  at  the  oxide-silicon  interface  corresponding  to  relaxation  times  (^)  of  100 
hours,  1  hour,  and  1  minute  respectively.  The  normal  component  of  the  surface  traction  is 
plotted  for  every  time  step,  just  like  the  outline  of  the  oxide  in  Fig.  1.  (Because  the  oxide 
shape  is  almost  identical  for  all  relaxation  times,  only  one  is  shown  in  Fig.  1.)  In  all  the  3 
stress  plots,  there  are  two  peaks  in  the  compressible  stress  region  (negative  value  range). 
The  early  peak  occurs  at  the  edge  of  the  nitride  mask  ( x  =  0.96//m).  This  peak  is  due  to 
the  highly  nonuniform  oxidation  rate  in  that  region.  As  time  progresses,  the  peak  shifts 
to  the  left,  further  into  the  nitride  mask,  and  gives  rise  to  a  late  peak.  As  expected,  stress 
decreases  as  the  viscosity  gets  lower. 

In  the  second  part,  we  repeat  the  same  simulations  with  the  silicon  nitride  layer  modeled 
as  an  elastic  material.  The  Young’s  modulus  and  the  Poisson’s  ratio  of  silicon  nitride  are 
assumed  to  be  3.29  x  1011  dynes-cm-2  and  0.266  respectively  (6} .  The  thickness  of  the 
nitride  layer  is  0.1/im.  The  final  shapes  of  the  oxide  and  the  nitride  layer  are  shown  in 
Fig.  3a,  3b,  and  3c.  As  we  can  see,  the  nitride  layer  bends  less  as  the  oxide  becomes  less 
viscous  and  flows  more  readily.  The  corresponding  values  for  peak  stress  are  2  x  10n, 
1.4  X  1011,  and  4.6  x  1010  dynes-cm-2  respectively.  The  first  two  values  are  unrealistically 
high  because  they  are  the  same  order  of  magnitude  as  the  elastic  moduli  of  oxide.  We  find 
that  in  general  we  need  those  stress  values  in  order  to  bend  the  nitride  mask  to  a  degree 
comparable  to  what  is  found  in  experiments.  To  keep  stress  down  to  a  realistic  level,  silicon 
nitride  must  deform  elastoplastieally  or  viscoelastically.  Unfortunately  we  don’t  have  data 
on  those  behaviors. 

4  Summary 

A  boundary  element  technique  has  been  developed  to  model  thermal  oxidation  of  silicon 
in  two  dimensions.  It  can  handle  a  wide  range  of  relaxtion  times  for  the  viscoelastic  flow 
of  oxide.  Simulations  of  some  simple  structures  have  been  demonstrated.  Preliminary 
results  indicate  that  it  is  inadequate  to  treat  silicon  nitride  as  an  elastic  material.  A 
comprehensive  characterization  of  silicon  nitride  is  clearly  needed. 


where  E  is  the  Young’s  modulus,  G  the  shear  modulus,  and  u  the  Poisson’s  ratio.  The 
fundamental  solutions  are  functions  of  time;  they  are  evaluated  at  the  end  of  a  time  step 
of  size  At.  Ka  and  Kp  decay  exponentially  with  a  time  constant  of  ra  and  Tp  respectively 
when  the  load  is  removed. 
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Abstract 

In  this  paper,  the  presence  of  delay  in  a  job 
shop  is  addressed.  We  show  that  delay  is  an 
important  consideration  in  many  manufactur¬ 
ing  systems  that  are  modeled  as  continuous  flow 
processes.  A  scheduling  policy  for  a  job  shop 
with  delays  is  then  derived  using  theoretical  ar¬ 
guments  and  heuristics. 

1  INTRODUCTION 

It  is  well  known  that  the  optimal  solution  of 
the  job  shop  scheduling  problem  is,  in  general, 
NP-hard  [2],  Except  for  a  few  problems  un¬ 
der  very  specific  conditions,  no  computationally 
tractable  solution  for  optimization  can  be  found. 
Due  to  this  formidable  computational  complex¬ 
ity,  which  necessitates  the  use  of  static,  over¬ 
simplified  models,  traditional  job  shop  schedul¬ 
ing  approaches  have  not  proven  satisfactory  in 
practice. 

The  approach  proposed  in  [l],  which  in  turn 
is  a  natural  extension  of  [3],  makes  use  of  a  hier¬ 
archical  control  structure  to  remedy  these  prob¬ 
lems.  A  high  level  controller,  similar  to  what 
described  in  [3],  works  at  long  time  scales  and 
deals  only  with  work  stations  (work  centers).  It 
treats  the  production  process  as  a  continuous 
material  flow.  Its  objective  is  to  control  the  flow 
over  a  long  time  horizon  so  that  the  demand  is 
satisfied  as  closely  as  possible  and  inventories 
are  kept  low,  while  keeping  the  system  within 
production  rate  capacity  constraints. 


The  actual  loading  of  individual  parts  into 
machines  is  left  to  lew  level  controllers  which 
work  at  shorter  time  scale.  The  low  level  deals 
only  with  single  work  stations  which  have  far 
fewer  machines  than  the  whole  job  shop.  The 
low  level  attempts  to  fulfill  the  production  goal 
determined  by  the  high  level  controller.  In 
this  way,  the  two  level  controller  can  avoid  the 
formidable  computation  requirements  encoun¬ 
tered  in  traditional  approaches.  Further,  it  dy¬ 
namically  adjusts  the  production_to-Cope  with 
real-time  events. 

While  the  two-level,  continuous  flow  model 
does  simplify  the  job  shop  scheduling  problem, 
it  comes  with  a  hidden  cost,  namely  that  the  dif¬ 
ferential  equations  representing  the  system  must 
often  include  delay.  To  see  this,  notice  that  any 
work  station  that  typically  processes  many  parts 
at  a  time  (i.e.  where  the  number  of  total  parts  in 
processes  is  much  greater  than  l)  will  have  av¬ 
erage  interarrival  times  that  axe  much  less  than 
the  processing  time  for  a  single  part.  For  such  a 
system,  the  time  parts  spend  in  the  system  can¬ 
not  be  ignored  and  thus  delay  must  be  explicitly 
included  in  the  formulation. 

In  this  paper,  which  is  a  summary  of  the  work 
to  appear  in  [8],  we  analyze  the  high  level  con¬ 
troller  for  systems  with  delay.  In  Section  2  we 
look  at  some  examples  of  manufacturing  systems 
with  delay  and  show  that  a  rather  large  class 
of  manufacturing  systems  require  delay  formu¬ 
lations.  In  Section  3  we  then  show  that  a  delay 
system  can  be  approximated  by  a  system  of  first 
order  differential  equations  without  delay.  We 
use  the  results  of  [6]  and  let  our  approximation 
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Figure  1:  A  single  work  station  with  delay 


delay  =r  '-*OH  delay=0 


X2 


jO— 


get  arbitrarily  good  to  arrive  at  a  solution  for  the 
optimal  control.  Due  to  the  difficulty  of  com¬ 
puting  the  optimal  value  function,  we  next  ex¬ 
plore  a  suboptimal  strategy  based  on  quadratic 
approximations  to  the  value  function.  Finally, 
conclusions  are  presented  in  Section  4. 

2  The  Importance  of  Delay 
in  Manufacturing  Sys¬ 
tems 

We  mentioned  that  delay  arises  in  manufac¬ 
turing  systems  that  work  on  many  parts  at 
one  time.  We  will  now  examine  this  phe¬ 
nomenon  more  closely  and  also  try  to  indicate 
in  what  ways  delay  introduces  difficulties  into 
the  scheduling  problem. 

Firstly,  let  us  point  out  that  introducing  de¬ 
lay  does  not  necessarily  complicate  the  control 
problem.  Consider,  for  example,  a  single  work 
station  with  delay  as  shown  in  Fig.  1.  where  x(t) 
is  the  inventory  in  the  buffer,  r  is  the  delay  (pro¬ 
cessing  time),  u(t)  is  the  loading  rate,  which  is 
bounded  and  the  bound  itself  is  a  random  vari¬ 
able  ([3]),  d(t)  is  the  demand  rate,  which  is  as¬ 
sumed  to  be  deterministic  and  known.  The  dy¬ 
namics  of  this  system  can  be  modeled  as 

*W  =  u(f  -  r)  -  d[t)  (1) 

By  simply  defining  x[t)  =  x[t  +  r)  and  d(t)  = 
d{t  -1-  r),  both  of  which  can  be  determined  com¬ 
pletely  at  time  t,  we  can  see  this  problem  is  no 
different  than  the  non-delay  problem.  We  sim¬ 
ply  use  x[t  +  r),  rather  than  x(£),  as  the  current 
state  and  solve  the  problem  as  though  there  were 
no  delay. 


Figure  2:  A  two  stage  system 

Unfortunately,  delay  cannot  always  be  han¬ 
dled  so  simply.  For  example,  consider  the  simple 
two  stage  system  depicted  in  Fig.  2. 

The  system  is  described  by 


*i(0  =  (t-rx)-u2(t)  (2) 

i2(t)  =  xi2{t)-d[t)  (3) 

0  <  ii(t)  (4) 

0  <  ui(0<ai(i)  (5) 

0  <  “2  (t)  <  cr2(t)  (6) 


Suppose  the  constraints  for  uj  and  t/2  depend  on 
some  random  processes  (e.g.  the  machine  state). 
We  must  determine  and  u2  based  on  the 
present  constraints  yet  the  value  of  the  future 
inventory,  Xi(i),  depends  on  both  the  present 
ui (0  &nd  the  future  u2(t),  the  constraints  on 
which  we  do  not  know. 

Another  example  where  delay  makes  the  prob¬ 
lem  more  complex  is  in  the  scheduling  of  a  reen¬ 
trant  job  shop.  A  reentrant  job  shop  is  one 
where  parts  visit  the  same  work  station  several 
times  [7],  A  simple  reentrant  job  shop  is  shown 
in  Fig.  3. 

New  parts  are  processed  by  the  work  station 
then  go  back  to  the  same  work  station  for  a  sec¬ 
ond  process.  After  the  second  process  is  finished, 
they  leave  the  system.  There  are  buffers  after 
the  first  and  second  processes  whose  levels  are 
denoted  by  and  x2  respectively.  Suppose  the 
processing  time  for  the  second  process  is  negligi¬ 
ble.  We  then  get  the  same  system  equations,  (2) 
and  (3).  The  only  difference  now  is  that  the  con¬ 
straints  on  ui  and  u2  are  also  coupled,  namely, 
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Figure  3:  A  reentrant  job  shop 

aiui(£)  +  <*2“2(£)  <  a(t)  for  some  ai,a2  and 
a.  This  further  complicates  the  control.  Thus, 
a  single  reentrant  work  station  with  delay  also 
cannot  be  trivially  handled. 

In  the  next  section,  we  expand  on  the  ideas 
suggested  by  these  examples  and  define  the  con¬ 
trol  problem  in  exact  terms. 

3  Solution  For  Delay  Sys¬ 
tems 

To  demonstrate  our  solution  technique  more 
clearly,  we  first  investigate  the  simple  problem 
described  in  (2)  to  (6).  The  technique,  however, 
is  extendible  to  more  complex  systems. 

The  objective  functional  is 

.€T.)/s(ll,I!)*  (7) 

here  fi(a)  is  a  polyhedron  defined  by  (5)  and  (6) 
and  j(-)  is  some  function  of  Zj,  and  x3.  Without 
delay,  this  is  the  same  formulation  as  in  [3]. 

At  time  t,  the  parts  in  the  first  process  that 
were  loaded  between  t  -  rx  and  t  will  contribute 
to  the  future  inventory  and,  therefore,  should 
become  part  of  the  current  state.  Unfortunately, 
the  problem  then  becomes  an  infinite  dimen¬ 
sional  one. 

In  order  to  overcome  this  difficulty  we  approx¬ 
imate  the  past  uj  by  a  finite  dimensional  first 
order  system.  We  then  let  the  approximation 
become  better  and  better  so  that  it  approaches 
the  original  system. 


Figure  4:  Diagram  of  two  systems 

Let  us  first  define  new  variables  yx  (t)  to  ym(f) 
through  the  following  equations. 

—  yi(0  =  ui(t)-yi(t) 

=  Vi(t)  -ifc(t)  (8) 

=  !/m— 1  (t)  ~  t/m  (t) 

The  initial  conditions  are  set  to  zero  at  —  oo 
and  we  assume  that  u^-oo)  =  0.  Eq.  (8)  de¬ 
fines  a  cascade  of  m  first  order  systems  with 
time  constant  ~.  Its  input  and  output  are 
ui(t)  and  ym(f)  respectively.  As  a  motivation 
for  using  (8),  note  that  its  transfer  function 
is  1/(1  +  3r/m)m  which  yields  the  well  known 
limit  e~,T,  the  transfer  function  of  a  delay  r,  as 
m  —*  oo. 

Now  define 

Xl{t)  =  ym(t)~u2{t)  (9) 

x=(f)  =  U2{t)-d{t)  (10) 

Combine  (8)-(  10) ,  we  obtain  a  new  system. 
The  diagrams  of  this  system  or  the  original  sys¬ 
tem  defined  by  (2)  and  (3)  can  be  drawn  as  in 
Fig.  4. 

The  only  difference  between  the  system  de¬ 
fined  by  (2)  and  (3)  and  the  system  defined  by 
(8)-  (10)  is  the  first  box  Kt (s).  For  the  first 
system,  it  is  a  delay  element  with  delay  r.  For 
the  second  system,  it  is  a  linear  system  defined 
by  (3).  If  we  can  show  that  for  the  same  U|,  uj 
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Figure  5:  Compare  two  integrals 


and  d,  the  output  of  the  first  system  approaches 
the  output  of  the  second  one,  then  we  can  es¬ 
tablish  an  equivalence  between  the  two  systems. 
By  superposition,  it  is  sufficient  to  show  that 
the  integral  of  ym(£)  approaches  the  integral  of 
ui(t  —  r)  as  m  goes  to  infinity. 

To  show  this  result,  we  compare  the  two  sys¬ 
tems  shown  in  Fig.  5.  In  Fig.  5,  K i(s)  is  the 
system  defined  by  (8)  and  Jf2(s)  is  a  pure  delay 
of  r.  We  will  prove  that  Zj  (t)  approaches  z2(t) 
uniformly  in  £  as  m  — +  00.  First  we  need  the 
following  lemma  which  is  similar  to  [6]. 

Lemma  1  If  u(t)  is  differentiable  with  |u(£)|  < 
K  for  all  t  €  (—00,  +oo),  then 

lim  sup  ||si(£)-32(£)||  =  0 

00  tS(-oo,  +  oo) 


Proof:  ( see  [ 8 ]). 

The  integrals  in  Fig.  5  start  from  -00. 
Since  the  initial  conditions  are  zero  at  —  co  and 
u(— 00)  equals  zero,  we  can  switch  the  linear  op¬ 
erators  K\ ,  K2  with  the  integrators,  1/s.  Be¬ 
cause  u(t)  is  bounded,  the  integral  of  u(f)  has  a 
bounded  derivative.  Combining  these  facts  with 
Lemma  1,  we  obtain  the  following  lemma. 

Lemma  2  If  u(t)  £  L\  and  is  bounded,  then 

lim  sup  ||z1(£)-z2(£)||  =0 

—  oo.+oo) 


Using  Lemma  2 ,  we  see  that  the  output  of  the 
system  defined  by  (8)- (10)  approaches  the  out¬ 
put  of  the  system  defined  by  (2)  and  (3).  There¬ 
fore,  if  ux  is  optimal  for  the  first  system,  it  will 
also  be  optimal  for  the  second  one. 

We  will  now  consider  the  optimal  control 
for  the  system  defined  by  (8)- (10).  Define 
1  =  [xx  x2  yx  ...ymJ',  u  =  [ux  r2]'.  This 
system  can  be  written  in  a  compact  form  as 

x  =  Ax  +  Bu  +  Cd  (11) 

Using  the  same  approach  as  in  [3],  it  can  be 
shown  that  the  optimal  control  u  for  the  problem 
defined  by  (11)  and  the  constraints  (4)-(6)  can 
be  obtained  by  solving 

min  VxJ*(x,  a)Bu  (12) 

u60(a) 

Eq.  (12)  can  be  rewritten  as 

.  ,dJ*  dT,  dJ‘ 

mm  - - - — u2  +  — — ux 

uGn(a)  dx\  <5x2  dyi 

Where  J*[x,  a)  is  the  optimal  cost  to  go.  Unfor¬ 
tunately,  J *  remains  unknown  and  is,  even  for 
simple  problems,  difficult  to  compute.  There¬ 
fore,  we  seek  an  approximation  for  J *.  Ex¬ 
perience  ([4], [5])  shows  that  satisfactory  results 
can  be  obtained  with  relatively  crude  approx¬ 
imations  for  J* .  The  one  we  will  use  is  in  a 
quadratic  form  with  coefficients  that  are  func¬ 
tions  of  a,  namely, 


J'[x ,  a)  =  x'R[a)x  +  5(a)x 

(14) 

Then 

VxJ*(x,  a)  =  R[a)x  +  5(a) 

(15) 

Using  (15)  we  can  rewrite  (13)  as 

y 

min  [/?!  (x,  a)uj  +  #>(x,  a)u2] 

u60(a) 

(16) 

where 

m  2 

(x,  a)  =  Vi]  [<*)'J]  +  J2 

3= 1  j=i 

+  Pj{a) 

(17) 

Proof:  ( see  / 8 J). 


Letting  m  go  to  infinity  we  get 


fT  2 

P'  ~  J0  A(cri  a)tti(*-cr)icr+y^  fty(a)z, 

/=  1 

(18) 

Eq.  (16)  and  (18)  describe  the  sub-optimal  con¬ 
trol  law  for  our  system.  Note  that  instead  of 
a  very  complex  dynamic  program,  the  control 
(production  rates)  is  determined  by  a  linear  pro¬ 
gram  (16).  The  problem  is  further  simplified  due 
to  the  simple  structure  of  fl(a).  This  calculation 
can  easily  be  performed  in  real  time. 

The  terms  in  /?,•  are  easy  to  interpret.  Note 
that  Pi  is  a  function  of  x a  and  the  past  con¬ 
trol  u,.  The  second  and  third  terms  in  p{  are 
identical  to  the  terms  found  in  the  quadratic  ap¬ 
proximation  presented  in  [3]  for  a  non-delay  sys¬ 
tem.  Added  to  this  is  a  convolution  of  the  con¬ 
trol  ux(t)  with  some  weighting  function  /,(<r). 
If  a  constant  weighting  is  applied,  the  first  term 
would  be  the  integral  of  Ul(t)  from  t  -  r  to  t, 
which  is  the  number  of  the  parts  currently  un¬ 
der  processing.  Other  weighting  functions  are, 
of  course,  possible. 

Further  research  will  be  conducted  to  obtain 
the  correct  form  of  Pi.  Further,  since  a  real  man¬ 
ufacturing  system  is  much  more  complex  than 
the  system  described  here,  more  complex  sys¬ 
tem  models  will  be  investigated  which  will  take 
into  account  factors  such  as  reentrant  process 
and  time  varying  demand.  Finally,  numerical 
computations  and  simulation  experiments  will 
be  used  to  help  develop  solution  techniques. 


4  Conclusion 

In  this  paper,  we  examined  the  effects  of  delay 
in  manufacturing  systems,  we  proposed  a  gen¬ 
eral  model  for  a  network  of  work  stations  that 
includes  delay.  We  then  presented  a  technique 
for  analyzing  delay  systems  by  augmenting  the 
states  to  include  an  approximation  of  the  past 
control.  Using  quadratic  approximations  to  the 
optimal  value  function,  we  show  that  the  control 
takes  a  particularly  simple  form. 
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ments  (e.g.,  tolerance  better  than 
0.005"),  and  are  adaptable  (programmable) 
so  that  they  can  handle  different  types  of 
odd  components. 

As  a  result,  different  operations 
'insertions  of  odd  component  types)  can  be 
performed  by  different  robots,  but  the 
amount  of  time  required  for  a  given  opera¬ 
tion  depends  on  the  speed  of  the  robot 
that  performs  it.  Also,  each  robot  has 
different  configurations  (e.g.,  tools)  and 
inherent  capabilities  (e.g.,  accuracy  and 
repeatability) ,  which  results  in  different 
subsets  of  operations  that  each  robot  can 
handle  (with  nonempty  intersections  among 
those  subsets) . 

As  a  consequence,  not  only  does  the 
input  rate  of  part  types  into  the  system 
ha/e  to  be  determined,  but  also  the  deci¬ 
sion  of  where  to  send  each  part  for  each 
operation  (among  the  possible  alternative 
robots)  has  to  be  made. 

Such  systems  are  usually  justified 
economically  only  if  the  production  volume 
is  quite  high  (e.g.,  hundreds  of  thousands 
of  components  inserted  per  year)  and  the 
variety  is  high.  Because  of  their  flexi¬ 
bility,  they  art  ex  lected  to  meet  demands 
that  vary  in  tha  short  term  and  that  re¬ 
quire  high  utilizition.  The  work  presen¬ 
ted  here  aims  to  improve  system  perfor¬ 
mance  (e.g.,  to  lead  to  higher  throughput 
and  reduced  WIP  while  meeting  production 
demands) . 

Another  type  of  manufacturing  system 
is  comprised  of  conventional  and  advanced 
machining  centers.  The  latter  are  capar.lo 
of  performing  different  operations,  with 
varied  capabilities.  For  example,  some 
machining  centers  can  do  drilling  and 
milling  operations  that  otherwise  reouire 
two  different  conventional  machines.  Also 
there  are  3-  and  5-axis  machining  centers. 
The  latter  can  do  more  operations  than  the 
former  without  changing  the  part  fixtu- 
ring. 

As  in  an  electronic  insertion  system, 
the  scheduling  and  routing  problem  in  a 
system  of  several  machining  centers  is  not 
only  to  decide  on  the  input  flow  rate  of 
each  part  type,  but  also  where  each  opera¬ 
tion  should  be  done  among  alternative 
machines  with  different  capabilities. 

Literature  Survey 

In  this  paper  we  present  a  method 
that  considers,  at  the  same  time,  two 
functions  —  short-term  scheduling  and 
routing  —  based  on  a  global  view  of  the 
system.  Many  references  consider  just  one 


of  these  functions.  For  exam'  'e,  Whitt 
(1986)  presents  a  method  whiun  can  be  used 
just  for  the  local  routing  decisions. 
Although  him  paper  develops  generic 
queueing  methodology,  we  use  his  results 
to  show  an  example  of  local  routing  consi¬ 
derations. 

By  local  rcuting  decisions  we  refer 
to  a  situation  by  which  a  customer  (or  a 
part)  has  to  join  one  of  several  queues. 
These  queues  represent,  for  example,  the 
input  buffers  to  workstations.  The  alter¬ 
native  queues  are  those  of  the  alternative 
workstations  that  can  perform  the  next 
operation  on  a  part,  which  has  just  fi¬ 
nished  a  particular  operation. 

Whitt  shows  that  in  some  cases,  the 
system  average  delay  is  not  always  mini¬ 
mized  by  customers  joining  the  queue  that 
minimizes  their  own  individual  expected 
delay.  This  result  suggests  that  deci¬ 
sions  should  be  made  only  when  taking  a 
global  view  of  the  system. 

Routing  is  treated  in  papers  by  Hahne 
(1981),  Tsitsiklis*  (1981),  and  Seidmann 
and  Schweitzer  (1984) .  Hahne  and  Tsitsik- 
lis  deal  with  only  two  choices  and  ma¬ 
chines  whose  randomness  is  due  to  failure 
and  repair.  Seidmann  and  Schweitzer  have 
many  choices,  but  the  randomness  is  due  to 
variations  in  processing  times.  In  all 
cases,  the  full  system  is  not  considered. 
Instead,  o..iy  one  decision  point  is  con¬ 
sidered,  and  decisions  are  made  on  a  pure¬ 
ly  local  basis. 

By  contrast,  we  consider  the  whole 
system  and  do  not  treat  local  conditions 
in  detail.  This  suggests  that  a  hierar¬ 
chical  decision  policy,  in  which  both 
kinds  of  decisions  —  local  and  global  — 
are  made  separately,  may  be  appropriate. 
The  local  decisions  should  be  made  in  a 
way  that  is  consistent  with  the  decisions 
made  on  a  global  basis. 

Outline  of  Paper 

Section  2  states  the  problem.  Sec¬ 
tion  3  contains  our  solution,  which  is 
based  on  dynamic  programming.  Section  4 
describes  some  numerical  examples  and 
simulation  results.  Conclusions  and  new 
research  directions  are  discussed  in  Sec¬ 
tion  5. 

2.  PROBLEM  STATEMENT 

Section  1  describes  two  situations  in 
which  short-term  scheduling  and  routing 
decisions  are  required.  In  this  section 
we  represent  such  manufacturing  systems 

with  a  mathematical  model. 


1 


£  am  for  every  machine  m.  (5) 


3.  SOLUTION 


y  yk  =  V  y*n 
5T  nm  J  run 

and  all  part 


for  all  k  *■  Kn 
types  n, 


(6) 


where  Kn  is  the  name  of  the  first  ope¬ 
ration  performed  on  parts  of  type  n.  De¬ 
note  by  fi  (at)  the  set  of  all  y  flow  rates 
that  satisfy  (4)  -  (6) . 


Following  the  usual  dynamic  program¬ 
ming  practice,  define 

J(x,  a,  t)  » 


min  E  I 

ycQ(a) 


J  g(x(s) )ds  x(t)=x,a(t)=a 


(8) 


Note  that  0(a)  is  a  random  set.  As 
machines  fail  and  are  repaired  the  instan¬ 
taneous  capacity  changes.  The  rates  that 
material  flows  into  the  system  must  change 
as  0(a)  changes,  as  well  as  in  anticipa¬ 
tion  of  these  changes. 


This  function  satisfies  the  Bellman  equa 
tion  (Bertsekas,  1976) ,  which  takes  the 
following  form: 

0  =  min  (  g[x(t)  ]  +  Z  y''n-dn] 

yefilal  1  ^ L  n3xnl  m  "*  nJ 


Cost  Function 

We  seek  a  policy  that  minimizes  a 
cost  of  the  form 

J(x0,  a0,  0)  = 

E  [  J  g(x(s)  )ds  |  x(0)«xa,a(0)=aa  j  (7) 

in  which  T  is  the  short  term  period,  such 
as  an  eight  hour  shift  and  g ( • )  is  a 
positive  convex  function.  We  assume  the 
cost  function  does  not  reflect  true  costs, 
but  instead  is  chosen  to  lead  to  desirable 
behavior.  Thus,  the  details  of  g(-) 
are  not  important.  In  Section  3  we  de¬ 
scribe  an  approximation  method  which  uses 
only  certain  features  of  the  cost  func¬ 
tion. 

Dynamic  Programming  Formulation 

The  optimization  problem  can  be  writ¬ 
ten: 

minimize  J(xQ,  a0,  0) 


+  |f  +  P,  t)  }  .  (9) 

This  equation  has  the  following  in¬ 
terpretation:  we  seek  a  function  J(x,a,t) 
such  that  the  values  of  y(x,a,t)  £  0(a(t)) 
that  minimize  the  right  hand  side  of  (9) 
cause  that  expression  to  be  zero.  This  is 
a  nonlinear  partial  differential  equation 
which  we  cannot  expect  to  have  an  analytic 
solution.  (However,  in  the  case  of  a 
single  part  type  and  a  single  machine, 
Akella  and  Kumar  (1986)  were  able  to  find 
a  closed  form  solution.) 

If  (9)  has  a  solution,  the  optimal 
control  y  satisfies  the  following  linear 
programming  problem.  Note  that  the  cost 
coefficients  are  time-varying  . 

rain  f  I  y<B] 

n3xn  lii  “J 

subject  to 

yen(«) 


subject  to  dynamics  given  by  (1)  and  (2) 
and  y  £  0(a). 

Comparison  with  Kimemia  and  Gershwin 

Kimemia  and  Gershwin  (1933)  formu¬ 
lated  an  optimization  problem  in  terms  of 
u  of  equation  (3) .  This  formulation  is 
correct  when  there  are  no  route  choices 
except  among  identical  machines.  However, 
they  assumed  that  they  could  ignore  (6) 
even  when  route  choice  existed,  and  then 
determine  y  from  u  after  solving  the  prob¬ 
lem.  This  assumption  is  not  correct;  the 
above  formulation  is.  Without  (6) ,  the 
choice  of  routes  achieved  n.ay  not  be  feas¬ 
ible,  and  (3)  would  not  necessarily  hold. 


It  is  important  to  recognize  that 
this  is  a  feedback  control  law  since  J  and 
n  are  functions  of  x  and  a.  The  solu¬ 
tion  y  is  therefore  a  function  of  x  and 
a  . 

Note  that  J  is  positive  since  it  is 
the  expected  value  of  the  integral  of  g,  a 
positive  quantity.  Note  also  that  feed¬ 
back  law  (10)  minimizes 


dJ  _  3J  T  y 

n  m 


(ID 


while  a  is  constant.  This  is  because  y 
appears  in  (9)  only  in  the  same  term  in 
which  it  appears  in  (11).  If  ot  remains 
constant  long  enough,  and  there  is  a  y  £ 
0(a))  such  that:  (11)  is  negative,  then  J 


JttlNYfcWSfcVV VV\YW%WU*VV  A  WL>VV  wN  1VLV -V.VA  A  WV  A.'-  A  A  A  A  A  L«s  A  A  AAA  AAA  V*  A  i>  A 


cated  is  which  of  the  following  conditions 
that  determine  the  values. 


>  0  (Regions  I  and  III) 
=»  v1  =  0 ,  v1  =0 

Ml  '  1  13 

>  0  (Regions  I  and  II) 


(A) 


3J 

3x, 


3J 


*  4  =  °'  4  "  0  (B) 

<  0  (Regions  II,  IV,  V,  and  VI) 


(C) 


<  0  (Regions  III,  IV,  V,  and  VI) 


Y  ,  *  Zf 


(D) 


3J 


g“-  <  0  and  >  o  (Region  II) 


'13 


3x2 

t!3 


> 

3x, 


0  and  < 
3xl 

^  y2  -  % 

23  T23 


0  (Region  III) 


(E) 


(F) 


If  both  derivatives  are  negative 


(Regions  IV  and 
are  already 


‘22 


v)  ,  y|,  and 

determined. 

1  fi  -  1,  2)  , 


The  remaining  variables,  yj3  (i  * 
minimize 

M.vi  +  iJv2 
3x,Yl3  +  3x2y23 

subject  to  (16) .  The  solution  is 


3J 

3J 

3x2 

T!3 

3X| 

2 

«3 

=> 

J23 

3J 

JL 

3J 

3x2  ~ 

T‘l3 

3x, 

i 

«3 

=»  y1 

*  13 

3 

y2 

‘23 


0. 


(21) 


(G) 


(H) 


In  each  of  these  regions,  the  control 
V1  moves  the  state  x„  through  the  dyna- 

ran 

mics  [(18)  and  (19)].  The  state  moves  to 
a  boundary  and  then  to  another  region. 
However,  there  is  one  exception.  In  both 
Regions  IV  and  V  the  state  moves  toward 
the  common  boundary,  which  is  given  by 

{  +  ft  '  t  =  0  (Region  VI),(I) 


If  we  follow  rules  (G)  and  (H) ,  the 
state  will  move  back  and  forth  across  the 
boundary  in  an  unrealistical  manner.  This 
is  called  chattering ,  It  occurs 
because  the  problem  is  singular ,  and 
a  remedy  is  suggested  by  Gershwin,  Akella, 


and  Choong  (1985).  A  strategy  is  found 
which,  when  x  reaches  Region  VI,  keeps  x 
in  Region  VI.  That  is,  it  maintains  (I). 
It  does  this  by  determining  y* 

nm 

which  minimizes  (21)  subject  to  (16)  and 


A _ /.L.  J£ 1  i£\  =  0 

dt(Tj33x2  T>  3x,  j  °* 


(22) 


This  is  simplified  by  assuming  that  J 
is  quadratic: 


J  =  ^xrA(a)x  +  b(a)Tx  +  c(a) 
Then 


i£ 

3x, 


A,,x,  +  A.,x,  +  b. 


and 


=  A21X2  +  A22X2  +  b2‘ 


If 


3J  ,  n 
5x7  0 


and  y  is  chosen  so  that 

r  i  33  i  jj  ]  E  n 


(23) 


(24) 


(25) 


(26) 


then 


(A2,x,  +  A^Xj  +  b2) 

T23 

-  +  (A„x,  +  A,2x2  +  b,)  =  0.  (27) 

‘  13 

Since  this  is  true  for  more  than  just 
one  instant,  its  first  derivative  with 
respect  to  t  is  also  0.  That  is, 


t2  { AzX|  + 
T23 

A22x 

2  +  b2] 

_  _L 

t!3 

(A„X 

i  +  AI2x2  + 

b.)  =  0, 

(28) 

or, 

f 

# 

A2.K+y;3-d0+A4 

YL+y23*dO 

+  b2J 

i 

T'i3 

(AnC^Yla^: 

!+Au( 

y22+y23-d:' 

)+b,^=°, 

(29) 

7 


From  (C)  and  (D) , 

-H  Aa(  oe./T'+y^-d,  )+A22(  «j/ii+34-d2  )+b2 ) 

723  ^ 

-+f A11[«1/T!1+y'13-d1)+A12(a2/i2+y223-d2)+b1]=o.  (3  0) 
M3*" 

Now  (16)  (as  an  equality)  and  (30) 
are  two  equations  in  two  unknowns,  y13  and 

y*  .  The  solution  is 


After  x  arrives  at  Region  VI,  it 
stays  in  Region  VI  if  yjf  and  y^  are 
given  by  (C)  and  (D)  and  yj3  and  y^  are 

given  by  (31)  and  (32).  Chattering  is 
avoided. 

5.  CONCLUSIONS 

This  paper  presents  an  extension  to 
the  earlier  Kimemia  and  Gershwin  work  to 
add  a  real-time  routing  calculation  to 
real-time  scheduling.  Thus  this  model  can 
be  used  for  many  more  types  of  manufactu¬ 
ring  systems. 

Future  work  will  include  the  develop¬ 
ment  of  local  operational  rules  which 
follow  the  system  routing  decisions  calcu 
lated  here,  and  extensive  simulation  of 
various  types  of  industries  to  further 
demonstrate  the  use  of  this  work. 
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